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The coordinate space formulation of the Hartree-Fock-Bogoliubov (HFB) method enables self- 
consistent treatment of mean-field and pairing in weakly bound systems whose properties are af- 
fected by the particle continuum space. Of particular interest are neutron-rich, deformed drip-line 
nuclei which can exhibit novel properties associated with neutron skin. To describe such systems 
theoretically, we developed an accurate 2D lattice Skyrme-HFB solver hfb-AX based on B-splines. 
Compared to previous implementations, we made a number of improvements aimed at boosting the 
solver's performance. These include: explicit imposition of axiality and space inversion, use of the 
modified Broyden's method to solve self-consistent equations, and a partial parallelization of the 
code. hfb-AX has been benchmarked against other HFB solvers, both spherical and deformed, and 
the accuracy of the B-spline expansion was tested by employing the multiresolution wavelet method. 
Illustrative calculations are carried out for stable and weakly bound nuclei at spherical and very 
deformed shapes, including constrained fission pathways. In addition to providing new physics in- 
sights, hfb-AX can serve as a useful tool to assess the reliability and applicability of coordinate-space 
and configuration-space HFB solvers, both existing and in development. 

PACS numbers: 21.60.Jz, 21.10.-k, 02.70.-c, 31.15.cj 



I. INTRODUCTION 

The important new focus in theoretical nuclear struc- 
ture research is to develop a coherent theoretical frame- 
work aiming at the microscopic description of nuclear 
many-body systems and capable of extrapolating into un- 
known regions. An important component in the theoret- 
ical landscape, and a crucial part of the theory roadmap 
Q,Q{, is the nuclear Density Functional Theory (DFT) in 
the formulation of Kohn and Sham [3(. Since the major- 
ity of nuclei in their ground states are superconductors, 
pairing correlations have to be taken into account. The 
resulting HFB or Bogoliubov de-Gennes equations can 
be viewed as the generalized Kohn-Sham equations of 
the standard DFT. The main ingredient of the nuclear 
DFT 4, 5] is the energy density functional that depends 
on proton and neutron densities and currents, as well as 
pairing densities representing correlated nucleonic pairs 

a- 

The unique structural factor that determines many 
properties of weakly bound nuclei is the closeness of the 
particle continuum. While the nuclear densities of bound 
nuclei eventually vanish at large distances, the wave func- 
tions of positive-energy states do not decay outside the 
nuclear volume, and this can be a source of significant 
theoretical difficulties. This problem is naturally over- 
come in the HFB method with a realistic pairing interac- 
tion in which the coupling of bound states to the particle 
continuum is correctly taken into account. In this con- 



text, particularly useful is the coordinate-space formula- 
tion of HFB [3, HI] . Another advantage of this method is 
its ability to treat arbitrarily complex intrinsic shapes, 
including those seen in fission or fusion. 

The HFB equations of the nuclear DFT represent the 
self-consistent iterative convergence schemes. Their com- 
putational cost can become very expensive, especially 
when the size of the model space - largely determined 
by the self-consistent symmetries imposed - or the num- 
ber of nuclear configurations processed simultaneously is 
huge. The advent of teraflop computers makes such large- 
scale calculations feasible, but in order to better optimize 
unique resources, new-generation tools are called for. 

A number of coordinate-space approaches to the nu- 
clear DFT have been developed over the years, and their 
performance strongly depends on the size and symme- 
tries of the spatial mesh employed [3, Gil El]. The 
ev8 code solves the Hartree-Fock plus BCS equations 
for Skyrme-type functionals via a discretization of the 
individual wave- functions on a 3D Cartesian mesh [12j |. 
assuming three symmetry planes. The ID hfbrad finite- 
difference code has been developed as a standard tool 
for HFB calculations of spherical nuclei [13[ . While lim- 
ited to the radial coordinate only, hfbrad allows very 
precise calculations, as the mesh step can be taken very 
low. The recently developed parallel 2D hfb-2d-lattice 
code [1J, [lji [lj| is based on B-splines; it can treat 
axially deformed nuclei, including those with reflection- 
asymmetric shapes. No coordinate-space 3D HFB nu- 



clear framework exists at present; however, a number of 
developments are under way, including a gene ral-purpose 
HFB solver based on wavelet technology [lj, LL8( ■ 

One traditional way of solving the HFB problem is 
to use the configuration-space technique, in which HFB 
eigenstates are expanded in a discrete basis, such as the 
harmonic oscillator (HO) basis [lQJ] . This method is very 
efficient and has been applied successfully in the large- 
scale calculations of nuclear properties |20J ■ However, the 
use of the HO basis is questionable in the limit of weak 
binding (because of incorrect asymptotic behavior of HO 
states) and at very large deformations, both requiring the 
use of unrealistically large configuration spaces. In both 
situations, a coordinate-state description is superior. 

The main objective of this paper is to develop a reli- 
able and accurate HFB axial solver, based on B-splincs, 
to study weakly bound nuclei and/or constrained en- 
ergy surfaces involving very large deformations. An HFB 
scheme based on such a concept has been implemented 
in hfb-2d-lattice by the Vanderbilt group [l4l. [TEL [l6| . 
An attractive feature of hfb-2d-lattice is that by tak- 
ing high-order B-splines one guarantees the correct rep- 
resentation of derivative operators on the spatial lattice 
0| . Unfortunately, the performance of hfb-2d-lattice 
is strongly CPU-limited. To speed up the calculations, 
and at the same time to improve the accuracy, we wrote 
a new 2D HFB B-spline code, called hfb-AX, in which 
a number of new features were introduced. Firstly, we 
incorporated space inversion as a self-consistent symme- 
try. Since reflection-asymmetric ground-state deforma- 
tions are present only in a handful of nuclei, this is not 
a serious limitation. Furthermore, we improved the it- 
erative algorithm by means of the modified Broyden's 
method. This resulted in significant convergence acceler- 
ation. We also made a number of smaller optimizations to 
the HFB solver which are described in the text. The per- 
formance and accuracy of hfb-AX has been very carefully 
tested against other codes. In short, we developed a fast, 
deformed coordinate-space HFB solver that can be used 
to carry out large-scale calculations on leadership-class 
computers, and it can also be invaluable when bench- 
marking the next-generation nuclear DFT tools. 

This paper is organized as follows. Section UT1 shortly 
summarizes the coordinate-space HFB approach in the 
cylindrical coordinate system and describes the numer- 
ical method used. In Sec. IIII1 numerical tests are pre- 
sented, with an emphasis on benchmarking against other 
existing HFB solvers. The examples include: (i) a two- 
center potential problem and a comparison with the mul- 
tiresolution wavelet method and the HO expansion tech- 
nique; (ii) comparison with the spherical finite-difference 
solver hfbrad for stable 120 Sn and drip-line Ni iso- 
topes; (iii) study of neutron-rich deformed 102 > 110 Zr and 
comparison with the axial solvers hfbtho and hfb-2d- 
LATTICE; and (iv) fission path calculations for 240 Pu. Fi- 
nally, Sec. HVI contains the main conclusions of this paper. 



II. THEORETICAL FRAMEWORK AND 
NUMERICAL METHOD 

A. HFB equation in cylindrical coordinate space 

The HFB equations in the coordinate-space represen- 
tation can be written as 0, [g, UM '■ 



dv r 



h(ra,r'a')-\ h(ra,r'a') 
h(ra,r'a') -h(ra,r'a') + A 



V> (1) (rV) 
</,(2)( r V) 



E 



V> (1) (rcr) 
V> (2) (ro-) 



(1) 



where (r, a) are the particle spatial and spin coordi- 
nates, h(ra, r'cr') and h(ra, rV) are the particle-hole (p- 
h) and particle-particle (p-p) components of the single- 
quasiparticle Hamiltonian, respectively, ipn ( r &) an d 
ipn (rer) are the upper and lower components of the 
single-quasiparticle HFB wave function, and A is the 
chemical potential. The spectrum of quasiparticle en- 
ergies E is discrete for \E\ < —A and continuous for 
\E\ > —A. With the box boundary conditions, the con- 
tinuum is discretized. In practical calculations, the p-h 
channel is often modeled with the Skyrme energy den- 
sity functional, while a zero-range S pairing interaction 
is used in the p-p channel. This choice is motivated by 
the fact that zero-range interactions yield the local HFB 
equations in coordinate space which are easy to solve. 

In the axially symmetric geometry, the third compo- 
nent of the single-particle angular momentum, fi, is a 
good quantum number. The HFB wave function can thus 
be written as ^„(r, H,g) where r=((f>, p, z), q=±\ de- 
notes the cylindrical isospin coordinates, and n=zt-^ , ±|, 
.. The corresponding HFB wave function can be 



±i 



written as uA 



.",(1) 



*n(T,n,q) = (2) '" 



e j(n -^i 2 U/^,T) 
Following the notation of Ref . [14j , we introduce 



'2it 



(2) 



U. 



(1,2) 



nflq 



D 



(1,2) 



{Pi*) = i>nnq(p,zA), 
(1,2), 



nQq(P> z ) = ^nClqiPiZ'l), 



(3) 



where wave functions U and D denote the spin-up ((7=™) 
and spin-down (cr=— -k) spinor components. Since the 
time reversal symmetry is conserved, one can consider 
only positive-fi values. In terms of these wave functions, 



the particle density p q (r) and pairing density p q (r) can 
be written as: 



*(*) = EEA^Sh 






/5g( r ) 



^E Eti<(p^)i 2 

+ Pi 2 i(p^)| 2 ], 

-EEA(^)^;(^) 



(4) 






= --E Ei<(^)<(^) 

n=| -e„>o 

+ D%( P ,z)dX(p,z)}. 

In the above equations, the sums are limited by the quasi- 
particle energy cutoff £ max , defining the effective range 
for a zero-range pairing force, and the angular projection 
cutoff $7 max . This way of truncating the continuum space 
is different from that in the spherical hfbrad code [13] 
which employs a maximum j max cutoff on the total single- 
particle angular momentum. This implies that even with 
the same energy cutoff, hfb-AX and hfbrad have differ- 
ent pairing phase spaces. We shall return to this point 
later in Sec. IIII1 



a 2D box z max and p max . These boundary conditions are 
important for the construction of derivative operators. 

In a given (f2, q) block, the HFB Hamiltonian in Eq. (Q]) 
can be expressed through the mean fields h and h with 
specified spin projections [141 ]: 
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The local Skyrme p-h Hamiltonian h has the usual form 
01: 

h g (r, a, a') = - V • -^ V + U g - iB q ■ (V x a) (7) 
2m* 

where m* is the effective mass, U q is the central p-h 
mean-field potential including the Coulomb term for pro- 
tons, and the spin-orbit potential with 



B„ 



1 



W q (Vp(r) + Vp q (r)), 



(8) 



where W q is the spin-orbit coupling strength. 

The pairing Hamiltonian h corresponding to the zero- 
range density-dependent 5 interaction can be written as 



TABLE I: Boundary conditions of HFB wave functions at 
p=0 and z=0 in hfb-AX. The single-particle states are labeled 
by iV quantum numbers. 



Q — | = even 

7T = +1 



fi - | = odd 

7T = +1 



Q — | = even 
tt = -1 



n 



odd 



i 

2 

7T = -1 



P=0 

mr, 

dp 

DL= = 



|p=o 



dp 
U\ p=0 



| P =o = 




dp 
D\ p=0 



P =o = 




t/| p=0 = 



2=0 



at/ 

D\ z= o 



z =0 =0 





U\ z=0 



z=0=0 
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dz U " 
U\ z=0 



o = 
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9[7, 



; = = 



dz 
D\ z=0 =0 



In the reflection-symmetric version of the code HFB-AX 
discussed in this paper, we assumed the space inversion as 
a self-consistent symmetry. Consequently, the quasipar- 
ticle wave functions are eigenstates of the parity operator 
V: 



V^n,U,q{p, Z, 4>) = ir^n^qip, ~Z, 4> + 1"), 



(5) 



i.e., parity tt = ± is a good quantum number. The pres- 
ence of conserved parity implies specific boundary condi- 
tions at ,o=0 and z=0 (see Table H]). We also apply the 
box boundary conditions at the outer box boundaries; 
namely, the wave functions are put to zero at the edge of 



h q (r,<T,a') = V q p q (r)F(r)5 ac , 



(9) 



where Vq < is the pairing strength, and the pairing 
form factor -F(r) depends on the form of pairing Hamil- 
tonian [22j ]: 



F(r) = 



1- 



1- 



Po 
2 Po 



— volume pairing 

— surface pairing 

— mixed pairing 



(10) 



where po— 0.16 fm~ 3 . The volume pairing interaction acts 
primarily inside the nuclear volume while the surface 
pairing generates pairing fields peaked around or outside 
the nuclear surface. As discussed in Ref. [23|, different 
forms of .F(r) can result in notable differences of pairing 
fields in diffused drip-line nuclei. 



B. B-spline technique in hfb-AX 

The lattice representation of wave functions and the 
HFB Hamiltonian used in this work closely follows that 
of hfb-2d-lattice described in Ref. [l6J]. In hfb-AX, 
the wave functions are discretized on a 2D grid (r Q , zp) 
with the M-order B-splines: 



,,(1.2) 



Wktt*q(pci>*f>) 



J2Br(p a )Bf(z,)C^ 



y nSl 7r g 



(ii) 



where C lJ is the matrix of expansion coefficients; a = 
1, . . . , N p and j3 — 1, . . . , N z . The four components of 
HFB wave functions ([2|) are thus represented in a ma- 
trix form. The derivative operators are constructed us- 
ing the Galerkin method. In the B-spline representation, 
the HFB Hamiltonian acts on 2D wave functions like a 
tensor, i.e., 



associated wave functions (\P*| and |^) are mapped into 
the bra and ket vectors, respectively: 



m s \H$\* afi )=E. 



(13) 



h"s*l>(p a ,Zf3) =ip'(py,zs). 



(12) 



The HFB equation is solved by mapping the Hamiltonian 
tensor into a matrix, which is then diagonalized. The 



J 



In the following, we give some details pertaining to the 
Hamiltonian mapping, because the mapping rule is dif- 
ferent with respect to the 16 individual blocks in Eq. ([6]). 
For the diagonal blocks, it is straightforward to map a 
tensor into a matrix (fc, k'): 



MTT) : 
Mil) : 

-h(TT) 

-h(ii) 



- a + N p N z 



k = (f3- l)N p + a 
k' = {5-l)Np + 1 

k = (f3-l)N p - 

k' = (S-l)N p + 7 + N p N z 

k = {/3- l)N p + a + 2 x N p N z 
k 1 = (S- l)N p + 7 + 2 x N p N z 

k=(/3- l)N p + a + 3 x N p N z 
k 1 = (5- l)Np + 7 + 3 x N p N z 



(14) 



Following the same rule, the bra and ket vectors are mapped into vectors with indexes k and k' , respectively. It 
is more complicated to map the off-diagonal blocks of (J6j> . For the four ket vectors, the mapped index k should 
correspond to different columns of the Hamiltonian blocks. For example, the index of the first ket component U^ 1 ' 
should correspond to the upper index of the first column of ([6]): 



([/(!)* £)(!)* J7(2)* £)(2)*) 




MTI) 

Mil) "A 

kn) 

Mil) 



MTT) 
MIT) 

-MTT) + a 

-MIT) 



\ 



MTI) 
MID 
-KU) 
-h(ii) + xj 



( UWjpgZf,) \ 

{/(2) 



(15) 



For the four bra vectors, the mapped index k' corresponds to different rows of the Hamiltonian. For example, the 
index of the first bra component t/W* corresponds to the lower index of the first row of ©: 



{\UW*{j^zs)\ £> (1) * C/ (2) * D<- 2 >) 



/ h lS (n)-x Mti) 
V 



^(4T) 
ft(TT) 

h(U) 



h(U)-X 

kn) 
ku) 



ths(U) MU) \ 



HIV MID 
-MtT) + A -h(U 

-MIT) -h(U) + Xj 



( u(1) \ 
D (l) 

C/(2) 
V D^ j 



(16) 



Following this scheme, the HFB Hamiltonian, represented by 16 tensor blocks, is mapped into a matrix of the rank 
4 x N p N z . The resulting HFB Hamiltonian matrix is diagonalized using lapack routines [24J. 



The Coulomb potential is obtained by directly inte- 
grating the Poisson equation, 



V 2 <Hp,z) = -47TC 2 p p (p,z), 



(17) 



where <fi is the Coulomb potential and p p is the proton 
density. The Poisson equation, discretized on a B-splinc 
grid, can be written in a matrix form. The boundary 
conditions at large distances are given by the multipole 
expansion of the Coulomb potential [16| . The gradient of 
the Coulomb potential at z=0 or p=0 is set to be zero be- 



cause of the symmetries imposed. These boundary con- 
ditions are incorporated in the lattice representation of 
the Laplace operator. 



C. Numerical speedup 

The size of the HFB Hamiltonian matrix depends on 
the box size R, the largest distance between neighbor- 
ing mesh points in the grid h (the B-spline grid is not 



uniform), and the order of B-splines M. Consequently, 
for large and refined grids, the diagonalization time be- 
comes a bottleneck. For example, for 120 Sn with a box of 
19.2 fm, a maximum grid size of h—0.6 fm, and M=13, 
the Hamiltonian matrix is about 0.45 GB of storage, and 
one diagonalization takes about 30 CPU minutes. 

Since the diagonalization of HFB matrices correspond- 
ing to different W blocks can be done independently, 
this part of hfb-AX can be parallelized using the Mes- 
sage Passing Interface. For 120 Sn with fi max =33/2 cut- 
off, about 70 processors are utilized. The precision of 
the derivative operators is crucial for the accuracy of 
the HFB eigenstates. By distributing diagonalization 
over many processors, one can perform calculations with 
larger boxes, denser grids, and higher-order B-splines. In 
addition, by taking advantage of the reflection symmetry, 
N z can be reduced; hence the rank of the Hamiltonian 
matrix is scaled down by a factor slightly less than 2. 
Consequently, the diagonalization process, which takes 
most of the execution time, can be speeded up by a fac- 
tor greater than 8. 

For the diagonalization, we employ the lapack dgeev 
routine [24j |. This routine diagonalizes a non-symmetric 
matrix using a QR algorithm. However, due to energy 
cutoff in HFB, it is not necessary to compute all eigen- 
vectors. For this reason, we modified dgeev so that it 
yields eigenvectors only within the required energy win- 
dow. By this way, the diagonalization time for 120 Sn is 
further reduced by one-third. 

To optimize the convergence of HFB iterations, we take 
the HF-BCS densities to warm-start the self-consistent 
process. Furthermore, instead of the commonly used 
linear mixing, we have implemented the modified Broy- 
den mixing, quasi-Newton algorithm to solve large sets 
of non-linear equations, for accelerating the convergence 
rate [25| • The iterative convergence can be estimated by 
the input and output difference at the m-th iteration, 



p( m ) 



T/( m ) 

v out 



V; 



(m) 



(18) 



where V is an ./V-dimensional vector containing un- 
knowns, and self-consistency condition requires that the 
solution V* be a fixed-point of the iteration: I(V*) = 
V* . For the linear mixing, the input at iteration m+1 is 
given as 



V; 



(m+l) _ y(ra) 



xF^ m \ 



(19) 



where a is a constant between and 1. In contrast, to 
estimate the next step, the modified Broyden mixing uti- 
lizes information obtained in previous Mb iterations. Re- 
cent implementations of this technique to the HFB prob- 
lem have demonstrated that the sufficient convergence 
can generally be obtained within 20-30 iterations [2y|. In 
HFB-AX the vector V consists of local densities and their 
derivatives at lattice mesh points: 

V = {Pq, T q , V • 3 q , p q , V%, V p p q , V z p q ) . (20) 

The dimension of V is thus 14x A p x N z . 



To demonstrate the performance of the method, we 
consider 22 in a 2D square box of i?=12fm, ft,=0.6fm 
and M=ll. The calculations were carried out on a Cray 
XT3 supercomputer at ORNL having 2.6 GHz AMD pro- 
cessors. Without reflection symmetry, one diagonaliza- 
tion takes about 15 minutes of CPU time. With reflection 
symmetry imposed, one diagonalization needs about 100 
seconds with the modified dgeev. The calculations with 
Broyden's mixing with Mb =7 are displayed in Fig. Q] 
The actual variation of binding energy is within 0.1 keV 
after 30 iterations. This is to be compared with linear 
mixing with a=0.6 which requires over 80 iterations to 
reach similar convergence. Both calculations show pre- 
cision limitations due to the numerical noise inherent to 
the mesh assumed. As in examples discussed in Ref. [26j , 
Broyden's method implemented in hfb-AX provides im- 
pressive performance improvements of HFB iterations. 
The numerical speedup is particularly helpful for heavy 
nuclei and for constrained calculations, which usually re- 
quire many self-consistent iterations. 
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FIG. 1: (Color online) Comparison between linear mixing 
(squares) and modified Broyden's method (circles) in hfb-AX 
for 22 0. The largest element of |p( m '| is shown as a function 
of the number of iterations m. See text for details. 



III. BENCHMARKING OF HFB-AX AND 
TYPICAL APPLICATIONS 



This section contains results of HFB-AX test results. 
First, the absolute accuracy of the one-body HF solver 
is tested using the adaptive multi-resolution method. 
Thereafter follows a series of calculations in which hfb- 
AX is benchmarked against the spherical coordinate space 
code hfbrad and the axial code hfbtho. In all those 
realistic tests, the Skyrme functional SLy4 [231 was used 
in the p-h channel, augmented by different density- 
dependent delta functionals ([5]) in the p-p channel. 



A. Two-center potential: comparison with HO and 
wavelet expansions 

The accuracy of the HFB-AX calculations in the p-h 
channel has been tested using an adaptive multiwavelet 
basis. To this end, we employed the madness framework 
|17| . The details regarding our particular realization of 
the wavelet basis expansion can be found in Refs. [1§,|28J. 

As a test case, we choose an axial two-center inverted- 
cosh potential: 

V(p, z) = V [f(p, z + Q + f( P , z-C)}, (21) 



where the inverted-cosh form factor is: 
f(p,z) ' 



1 + e - R «/ a cosh( x /p 2 + z 2 /a) 



(22) 



and Vb, Ro, and a are the potential's depth, radius, and 
diffuseness, respectively, and 2( is the distance between 
the two centers. A cross section of the potential used in 
our test calculations is shown in Fig. [2] along the z-axis 
at p=0. It is seen that the two potential's centers are 
well separated; hence the ability to predict a small par- 
ity splitting between the eigenstate doublets provides a 
stringent test for the B-spline Schrodinger equation solver 

Of HFB-AX. 




z(fm) 

FIG. 2: Fhe two-center inverted-cosh potential of Eq. (121[) 
as a function of z at p—0. The two centers are 15 fm apart 
(C=7.5fm) and F =-50MeV, 7? =2fm, and a=lfm. 



In addition to potential (|2ip we also considered the 
spin-orbit term in the usual Thomas form: 



Vso{p,z) = -iX 



2mc 



W(p,z)-(<rx V), (23) 



where Ao=5.0 and the numerical values of fundamen- 
tal constants were taken as h 2 /2m=20. 721246 fm 2 , 
fic=197.32696 MeV fm, and mc 2 =939. 56535 MeV. The 



corresponding one-body Schrodinger equation reads: 

h 2 



-^—V 2 + V(p,z) + V so (p,z) 
Am 



<p(p,z,<t>) =E(p(p,z,(f), 

(24) 
where ip(p, z, <f>) is a two-component spinor wave function. 



TABLE II: Ten lowest eigenvalues of the two-center potential 
(|21[) obtained using the HO, B-spline, and wavelet expansions. 
All energies are in MeV. For more details see text. 



State No. n* 



HO 

N sh =20 



HO 

N sh =30 



B-spline 
h=0.6 



Wavelets 



1 


1/2+ 


-22.23916 


-22.24008 


-22.24011 


-22.24011 


2 


1/2- 


-22.23816 


-22.23995 


-22.23998 


-22.23998 


3 


1/2+ 


-9.21514 


-9.22047 


-9.22050 


-9.22050 


4 


3/2" 


-9.20359 


-9.21256 


-9.21260 


-9.21260 


5 


1/2- 


-9.20359 


-9.21256 


-9.21260 


-9.21260 


6 


3/2+ 


-9.20589 


-9.21126 


-9.21129 


-9.21129 


7 


1/2+ 


-9.20589 


-9.21126 


-9.21129 


-9.21129 


8 


1/2- 


-9.19743 


-9.20590 


-9.20595 


-9.20595 


9 


1/2+ 


-1.70724 


-1.72402 


-1.72503 


-1.72514 


10 


1/2- 


-1.49218 


-1.52486 


-1.52672 


-1.52690 



Table [TT| displays the lowest eigenvalues of the two- 
center potential (|21[) obtained in three expansion meth- 
ods. In the HO expansion calculation, we took 
./V s /,,=20 and 30 shells of the spherical oscillator with 
?i(jJo=5.125MeV (as it turned out, the use of a stretched 
basis was not particularly advantageous). The size of the 
fT = i + Hamiltonian block is 121 and 256 for N sh =20 
and 30, respectively, i.e., the matrix size is more than 
doubled in the latter case. In the hfb-AX calculation 
with M=13, we used a square box of i?=25.2fm and 
/i=0.6fm. (The values of N sh =20 in HO and h=0.6im 
in HFB-AX are typical for realistic Skyrme-HFB calcu- 
lations.) In the wavelet variant [Ha, U^, the absolute 
accuracy was assumed to be 10" 5 . 

It is seen that the accuracy of B-spline expansion is ex- 
cellent, both for the absolute energies and for the parity 
splitting. The HO basis with N s h—20 performs rather 
poorly, especially for parity splitting and for the ener- 
gies of the highest (halo) states. That is, of course, 
to be expected for a two-center potential expanded in 
a one-center basis. It is only at iV s /j=30, not practical 
in large-scale DFT calculations, that a good agreement 
with wavelets and hfb-AX is obtained. In all variants, 
there is a perfect degeneracy oi P3/2 doublets with S!=i 
and £!=§ . 

The results with the inclusion of the spin-orbit term, 
which lifts the degeneracy between ^=5 and | levels, 
are given in Table IIIII It is seen that the general ex- 
cellent agreement between B-spline and wavelet variants 
holds, and that HO with N s h=20 performs rather poorly, 
especially for the halo state. Finally, the single-particle 
spectrum of a two-center potential is illustrated in Fig. [3] 
as a function of the inter-center distance. A transition to 
a dimer-like spectrum is clearly seen at distances greater 



TABLE III: Similar to Table [TT] except for the two-center 
double-cosh potential with the spin-orbit term. 



State No. fi" 



> 



HO 

N sh =20 



HO 

N sh =30 



B-spline 
h=0.6 



Wavelets 



1 


1/2+ 


-22.23916 


-22.24008 


-22.24011 


-22.24011 


2 


1/2- 


-22.23816 


-22.23995 


-22.23998 


-22.23998 


3 


1/2+ 


-9.43145 


-9.43659 


-9.43663 


-9.43662 


4 


3/2" 


-9.42314 


-9.43199 


-9.43203 


-9.43202 


5 


3/2+ 


-9.42561 


-9.43078 


-9.43081 


-9.43080 


6 


1/2- 


-9.41931 


-9.42783 


-9.42788 


-9.42788 


7 


1/2+ 


-8.77250 


-8.77825 


-8.77828 


-8.77828 


8 


1/2- 


-8.76475 


-8.77380 


-8.77384 


-8.77383 


9 


1/2+ 


-1.70727 


-1.72405 


-1.72506 


-1.72516 


10 


1/2- 


-1.49222 


-1.52490 


-1.52675 


-1.52693 




distance (fm) 

FIG. 3: (Color online) Eigenvalues of a two-center inverted- 
cosh potential with the spin-orbit term calculated in hfb-AX 
as functions of the distance between two centers 2(. 



MeV. For the pairing space, we adopted the commonly 
used equivalent energy cutoff of 60 MeV ;20j. As both 
codes are written in different geometry, the quasiparticle 
continuum is discretized differently in hfbrad and hfb- 
AX. In hfbrad, all partial waves with j<jmax = 33/2 
were considered, while in hfb-AX we imposed a cutoff on 
jz : ^max=33/2. For the sake of comparison, the pair- 
ing regularization option in hfbrad has been turned off. 
Also, we adopted the same fundamental constants as in 
Ref J2(i|: h 2 /2m = 20.73553 MeV and e 2 = 1.439978 MeV 
fm. 



TABLE IV: Results of spherical HFB+SLy4 calculations with 
volume pairing for 120 Sn using hfb-AX with different mesh 
size h and B-spline order M . The results of the precision 
radial code hfbrad are shown for comparison. The static 
proton pairing is zero. All energies are in MeV; h is in fm. 



E 



h=0.64 
M=ll 



h=0.6 
M=ll 



h=0.6 
M=13 



ft=0.15 

HFBRAD 



Etot 


-1018.304 


-1018.271 


-1018.356 


-1018.362 


E c 


347.636 


347.642 


347.577 


347.558 


Kin 


831.520 


831.512 


831.534 


831.520 


rpn 
^kin 


1339.516 


1339.553 


1339.562 


1339.598 


rpn 
f-Jpair 


-10.253 


-10.253 


-10.253 


-10.278 


rpn 
-&kin 


1329.263 


1329.300 


1329.309 


1329.320 


A n 


1.2451 


1.2451 


1.2451 


1.2450 


An 


-7.9950 


-7.9950 


-7.9950 


-7.9953 



Table IIVI displays various contributions to the binding 
energy E to t of 120 Sn, i.e., kinetic energy Ekin for protons 
and neutrons, Coulomb energy Eq, neutron pairing en- 
ergy Ep air , neutron pairing gap A„ and Fermi level A n , 
and the sum 



than 6 fm. Our tests nicely illustrate the ability of the B- 
spline technique (thus: hfb-Ax) to handle a two-center 
problem encountered, e.g., in fission or fusion [29j. 



B. Spherical limit: comparison with hfbrad and 

HFBTHO 

The performance of hfb-AX at the spherical limit can 
be assessed by comparing it against the accurate ID ra- 
dial coordinate code HFBRAD, based on a direct integra- 
tion of the system of coupled radial differential equations 
0; EDI ■ The tests have been carried out for the nucleus 
120 Sn, which is often used in benchmarking HFB solvers 
[la. |20| | in the limit of spherical shape and nonzero neu- 
tron pairing. 

The precision of HFB calculations in coordinate space 
is primarily determined by the size of the mesh used. We 
calculated 120 Sn with the fixed-box size (i?=19.2 fm) but 
with different mesh steps and B-spline orders. In our cal- 
culations, we took the volume delta pairing interaction 
with the pairing strength Vq=— 187.05 MeV fm 3 adjusted 
to the experimental neutron pairing gap of A„=1.245 



^hin -^A 



kin 



E 



pair 



(25) 



for neutrons. (As discussed in Refs. [30j, [3l|, [32J, while 
for zero-range pairing the individual values of Ekin and 
Epair are divergent with respect to the cutoff energy of 



the pairing window, their sum ([25)1 converges nicely and 
it is less sensitive to the actual treatment of discretized 
quasiparticle continuum.) The general agreement be- 
tween hfb-AX and hfbrad is excellent, in particular in 
the ft.=0.6fm, M=13 variant, with most quantities agree- 
ing within lOkeV. As expected, the largest difference is 
and E™„;„ due to a slightly different char- 



seen for E^ it 



pair 



acter of unbound spectrum in the two models; however, 
the sum Ekin is well reproduced by hfb-AX. 

The HFB results with mixed pairing obtained in hfb- 
AX, hfbrad, and HFBTHO are shown in Table [V] The 
pairing strength in hfb-AX was taken as Vq=— 284.29 
MeV fm 3 , as compared to -284.36 and -284.10 MeV 
fm 3 in hfbrad and HFBTHO, respectively [2(3]. In HF- 
BTHO calculations, 25 shells of transformed HO basis 
were used. Again, the agreement between hfb-AX and 
hfbrad is excellent, and the total binding energies ob- 
tained in the three methods agree within 12 keV. As ex- 
pected, the largest differences are seen for Ekin- In par- 



TABLE V: Similar to Table IIVI except for the mixed pair- 
ing interaction. hfb-AX results are compared with those of 
hfbrad and hfbtho of Ref. [13]. 



E 



ft=0.6 

M=13 



fo=0.1 

HFBRAD 



HFBTHO 



Etot 


-1018.795 


-1018.791 


-1018.777 


E C 


347.442 


347.400 


347.370 


EL 


830.856 


830.848 


830.735 




1340.675 


1340.668 


1340.458 


cm 


-12.491 


-12.467 


-12.467 


jpn 


1328.184 


1328.201 


1327.991 


A n 


1.2448 


1.2446 


1.2447 


An 


-8.0186 


-8.0181 


-8.0168 



ticular, hfbtho underestimates the neutron (proton) ki- 
netic energy by about 200 (100) kcV. This deviation is 
partly due to different representations of the kinetic en- 
ergy operator in the coordinate space and in the trans- 
formed oscillator basis, and partly due to the different 
continuum space (see the discussion above). 
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C. Weak binding regime: comparison with hfbrad 

Neutron-rich nuclei are unique laboratories of neu- 
tron pairing. In weakly bound nuclei, pairing fields 
are affected by the coupling to the continuum space, 
and this coupling can significantly modify pair distri- 
butions [3, [23, [33, [34|. In Sec. IIIIBI we demonstrated 
that HFB-AX performs very well for a stable spherical nu- 
cleus 120 Sn. To evaluate the performance of hfb-AX for 
weakly bound nuclei, in this section we discuss ground- 
state properties of even-even 84 > 86 > 88 > 90 Ni isotopes, which 
are expected to be weakly bound [10, Ha, [3a, [37| . 

In our test calculations, we adopted the surface pairing 
with strength adjusted to 120 Sn (Vo=-512.6MeVfm 3 ). 
Results of hfbrad and hfb-AX calculations are listed in 
Table I VII The nucleus 90 Ni is a drip-line system and 
its stability is strongly influenced by pairing. Indeed, 
it is only bound with surface pairing; with volume and 
mixed pairing, 90 Ni is calculated to have a positive neu- 
tron chemical potential (see Ref. [23j for more discussion 
concerning this point). The local particle and pairing 
densities of drip-line even-even Ni isotopes are shown in 
Fig. [4j A gradual increase of neutron skin as approach- 
ing 90 Ni is clearly seen. The proton density, on the other 
hand, is only weakly affected by the outermost neutrons. 

The systematic comparison between HFB-AX and HF- 
BRAD is given in Table fVTl For the binding energy, the 
agreement is very good, including the borderline system 
90 Ni. The neutron pairing energy increases as one ap- 
proaches the neutron drip line. This is consistent with 
the systematic behavior of pairing densities shown in 
Fig.H 

For 84 Ni, hfb-AX and hfbrad yield similar pairing 
properties and kinetic energies. However, with increas- 



FIG. 4: (Color online)Proton (top) and neutron pairing (bot- 
tom) densities calculated in hfb-ax for 84 > 86 ' 88 > 90 Ni. Proton 
pairing is zero. 



ing neutron number, the difference between the values of 
E% in obtained in the two models gradually grows, reach- 
ing over 1.6 MeV in 90 Ni. At the same time, the difference 
between E\ in values is smaller by an order of magnitude. 
The systematic difference between hfb-AX and hfbrad 
when approaching the neutron drip line can be traced 
back to different angular momentum truncations, i.e., 
pairing phase space structure. In hfbrad, the s.p. an- 
gular momentum cutoff is j m ax=33/2, while in the HFB- 
AX, the cutoff is done in terms of the s.p. angular mo- 
mentum projection and it is 5! max =33/2. Consequently, 
in HFB-AX contributions from high-j, low-fi continuum 
states, absent in hfbrad, are present. In Sec. IIII B[ we 
have shown that for the well-bound nucleus 120 Sn, this 
difference in the continuum phase space is insignificant. 
However, for nuclei close to the drip line, where the con- 
tribution from unbound states is far more important, the 
situation is very different. 

In order to quantify this point, we performed cal- 
culations for 90 Ni with j max =49/2 in hfbrad and 
f2 max =49/2 in hfb-ax. The results are displayed in Ta- 
ble IVIIl The variations in the proton kinetic energy be- 
tween various variants of calculations are small, suggest- 
ing that the kinetic energy operator is well represented 
by hfbrad and hfb-ax with the grids assumed. Also, 
the binding energy changes little, ~30keV, when the j- 
or J7-cutoff is increased. In the larger j-window, E^ in in 
hfbrad is reduced by about 1 MeV; the corresponding 
change in HFB-AX is much smaller, ~200 keV. This result 
indicates that the high-j continuum contributions play 
an important role in the structure of 90 Ni. A need for an 



TABLE VI: Comparison between hfb-AX and hfbrad with SLy4 p-h functional and surface pairing for drip-line nuclei 
84 ' 86,88,90 Ni. The same box 7?=19.2fm was used in both cases. The angular momentum cutoff was taken at j ma x=33/2 in 
hfbrad and f2 max =33/2 in hfb-AX. All energies are in MeV. 





84 Ni 




86 Ni 




88 Ni 




90- 


Nfi 


E 


HFBRAD 


HFB-AX 


HFBRAD 


HFB-AX 


HFBRAD 


HFB-AX 


HFBRAD 


HFB-AX 


Etot 


-654.919 


-654.899 


-656.915 


-656.955 


-658.215 


-658.193 


-658.877 


-658.856 


E c 


122.797 


122.806 


122.215 


122.228 


121.621 


121.640 


121.018 


121.056 


Ekin 


430.468 


430.460 


426.311 


426.330 


422.152 


422.206 


418.027 


418.204 


E?- 


1084.511 


1084.577 


1116.835 


1116.782 


1148.387 


1148.179 


1179.697 


1178.956 


jpn 
Impair 


-30.892 


-30.890 


-36.733 


-37.080 


-43.179 


-43.727 


-49.926 


-50.807 


jpn 


1053.619 


1053.687 


1080.102 


1079.702 


1105.208 


1104.452 


1129.771 


1128.149 


A„ 


1.485 


1.486 


1.613 


1.617 


1.742 


1.746 


1.862 


1.864 


An 


-1.455 


-1.454 


-1.062 


-1.068 


-0.709 


-0.718 


-0.399 


-0.417 



TABLE VII: Similar to Table |\T] except for 90 Ni with 
the angular momentum cutoff j ma x=49/2 in hfbrad and 
£} max =49/2 in hfb-AX. All energies are in MeV. 



HFBRAD 



HFB-AX 



Etot 

Ec 

E p 

kin 

f-^kin 
jpn 

jpn 
^kin 

A„ 

An 



-658.911 
121.038 
418.165 

1179.073 
-50.326 

1128.747 

1.860 

-0.410 



-658.881 
121.060 
418.233 
1178.843 
-50.892 
1127.951 
1.864 
-0.420 



appreciable angular momentum cutoff in the description 
of weakly bound nuclei, especially for surface-like pairing 
interactions, has been pointed out in Ref. [38j |. 



D. Deformed, weakly bound case: comparison with 

HFBTHO 

One of the objectives of hfb-ax is to precisely solve 
HFB equations for axially deformed nuclei, in particular 
at very large deformations and/or at the limit of weak 
binding. In this context, the neutron-rich Zr isotopes 
with A~110 are very useful testing grounds, as they are 
known to have very large prolate deformations [16j . In 
this section, we compare axial HFB-AX and hfbtho cal- 
culations for the nuclei 102 > 110 Zr which exhibit deformed 
neutron skin. 

Table IVIIII shows the results of deformed calculations 
for 102 Zr and 110 Zr with the same parameters as in spher- 
ical calculations for 120 Sn displayed in Table fVl In hf- 
btho, deformed wave functions were expanded in the 
space of 20 stretched HO shells. The binding energies 
in HFB-AX are greater by about 110-140 keV than those 
of HFBTHO. This is understandable as HFBTHO with 20 
shells also underestimates the binding energy of 120 Sn by 
about 150 keV [2(|. In Table I\TlIT1 we only show neutron 
pairing (proton pairing correlations in 110 Zr vanish due 
to a deformed proton subclosure at Z=40). It is grat- 



TABLE VIII: Results of deformed hfb-ax and hfbtho 
HFB+SLy4 calculations for 102 Zr and 110 Zr with mixed pair- 
ing. All energies are in MeV. The quadrupole moments are 
in fm 2 . 



-Zr 



'Zr 



HFB-AX 



HFBTHO 



HFB-AX 



HFBTHO 



Etot 


-859.649 


-859.540 


-893.983 


-893.840 


Ec 


231.149 


231.084 


226.758 


226.712 


^kin 


651.309 


651.099 


632.115 


631.882 


El- 
J -'kxn 


1202.050 


1201.990 


1368.206 


1368.201 


jpn 


-3.261 


-3.535 


-3.200 


-3.323 


jpn 


1198.789 


1198.455 


1365.006 


1364.878 


A n 


0.672 


0.700 


0.636 


0.652 


A-n 


-5.431 


-5.435 


-3.552 


-3.543 


Q20 


410.08 


411.31 


444.02 


443.90 


Q20 


638.19 


639.41 


788.32 


786.63 



ifying to see that energies and quadrupole moments of 
102 Zr and 110 Zr are very close in hfb-AX and hfbtho, 
in spite of fairly different computational strategies un- 
derlying these two codes. 



TABLE IX: Comparison between hfb-2d-lattice (sec- 
ond column) and HFB-AX (third column) for 102 Zr. Cal- 
culation parameters are the same as in Ref. [lj], i-e., 
Vo=-170MeVfm 3 , Q max =ll/2, #=12 fm, and N P =N Z =19. 
The results of hfb-AX with the same density functional 
but the standard box i?=19.2fm and larger pairing cutoff 
f2max=33/2 are displayed in the last column (hfb-ax'). All 
energies are in MeV. The mass rms radius R rm s is in fm. 





HFB-2D-LATTICE 


HFB-AX 


HFB-AX' 


Etot 


-859.61 


-859.19 


-859.25 


A 11 


-5.46 


-5.47 


-5.45 


Ap 


-12.10 


-12.05 


-12.06 


A„ 


0.31 


0.28 


0.43 


A p 


0.34 


0.37 


0.40 


■tirrns 


4.58 


4.58 


4.58 


fc 


0.431 


0.434 


0.435 



In order to compare HFB-AX with the Vanderbilt lattice 
code hfb-2d-lattice, we calculated the deformed nu- 
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cleus 102 Zr assuming the same parameters as in Ref. [14j . 
i.e., a fairly small box radius i?=12fm, coarse grid, and 
very low cutoff fi max =ll/2. The ground state of 102 Zr 
is reflection symmetric; thus, apart from the fact that 
the box of HFB-2d-lattice is twice as large as that of 
HFB-AX, the codes are supposed to produce the same re- 
sult. As seen in Table IIX1 this is almost the case: the 
difference for the binding energy is around 400 keV. We 
believe that this can be attributed to a slightly different 
structure of the discretized positive-energy continuum in 
the two codes. Indeed the average pairing gaps predicted 
in the two codes are ~30 keV apart. To confirm this, we 
performed another calculation for 102 Zr in a larger space, 
i.e., #=19.2 fm, h=0.6hn, and fi max =33/2. While the to- 
tal energy is only weakly affected, there is an appreciable 
increase in the pairing gaps. This result, together with 
the discussion of 90 Ni in Sec. lIIICl underlines the impor- 
tance of using large boxes and sizable pairing spaces for 
the description of neutron-rich systems. 
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FIG. 5: (Color online) Contour plots of proton (top) and 
neutron (bottom) density distributions in the (p, z)-plane for 
the deformed ground state of 110 Zr calculated in hfb-AX. The 
densities are in nucleons/fm -3 . 

Proton and neutron density distributions in 110 Zr are 
displayed in Figs. [5] (in two dimensions, to better show 



the deformed shape) and [6] (in logarithmic scale, to bet- 
ter show the asymptotic behavior). The appearance of 
the neutron skin beyond the nuclear surface is clearly 
seen. The density contours in Fig. O can be compared 
with the result of hfb-2d-lattice shown in Fig. 6 of 
Ref. [16j |. and there seems to appear a good agreement 
between the two sets of calculations. In particular, a 
small depression of density in the nuclear interior, due to 
shell effects, is present in both cases. Another interesting 
feature is the constancy of density diffuseness along the 
nuclear surface. The asymptotic behavior of nuclear den- 




P /z (fm) 

FIG. 6: (Color online) Particle (top) and pairing (bottom) 
ground-state densities in 110 Zr along p=0 and z—Q. The size 
of the box is #=19.2 frn. 

sities in Fig. [6] is consistent with general expectations [3], 
and the ratio p n (0, z)/p n {p — z, 0) is roughly constant at 
large distances. This indicates that densities are still well 
deformed in the region which is well beyond the nuclear 
surface. 



E. Large deformation limit: axial symmetric fission 
path of 240 Pu 

The advantage of coordinate-space calculations over 
HO expansion methods is apparent in the context of 
problems involving extreme deformations, which require 
the use of huge oscillator spaces or even a many-center 
HO basis. In this section, we study the axial, reflection- 
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symmetric fission path of 240 Pu, which has been investi- 
gated in many earlier works [4|. By carrying out precise 
HFB-AX calculations, one can assess the error on poten- 
tial energy surfaces, energies of fission isomers, and fission 
barriers obtained in commonly used HFB codes employ- 
ing an HO expansion technique. 

The HFB energy at a given mass quadrupole moment 
Q20 = (Q20) can be obtained by minimizing the routhian 
with a quadratic constraint [19j |: 



E' = E + C q ({Q 20 ) - Q 



211 J 



where 



?20 



2tt 



Ptot(p, z){2z 2 - p 2 )pdpdz, 



(26) 



(27) 



is the requested average value of the mass quadrupole 
moment and C q is the quadrupole stiffness constant. 
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FIG. 7: (Color online) Top: axial, reflection-symmetric fission 
path of 240 Pu (top) calculated with HFB-AX and hfbtho (in 
a spherical and stretched HO basis). Bottom: the difference 
between hfbtho and hfb-AX results (normalized to zero at 
the ground-state configuration). The minima and maxima 
of energy are marked: ground state, g.s.; first barrier, b(l); 
fission isomer, f.i.; second barrier, b(2). 

The constrained hfb-AX calculations for 240 Pu were 
performed in a box of i?=23.4fm and /i=0.65fm, using a 
B-spline basis with M=ll. (For a similar mesh size, the 
binding energy of 120 Sn in hfb-AX agrees with hfbrad 
within 50 keV.) The hfbtho calculations were carried 
out in a typical space of Ngh—20 spherical or stretched 



HO shells. The mixed pairing interaction is used with the 
pairing strength adjusted as in Sec. IIIIBl for 120 Sn. The 
results are displayed in Fig. [7] The spherical HO basis 
is unreliable for fission calculations, and the quality of 
hfbtho calculations with a stretched basis deteriorates 
gradually with deformation. That is, the energy error 
on the first barrier and fission isomer is ~100keV and 
^400 keV, respectively, and it reaches ~500keV inside 
the second barrier. These are significant corrections that 
can impact calculated half-lives for spontaneous fission. 



IV. CONCLUSIONS 

We developed a 2D coordinate space code hfb-AX us- 
ing the technique of basis splines. The high accuracy 
of the code has been demonstrated by benchmarking it 
against the multi-resolution wavelet method, the HO ba- 
sis expansion method, and the radial finite difference 
method. The tests have been carried out for spherical 
and very deformed configurations of stable and weakly 
bound nuclei. 

A significant numerical speedup of the code makes it 
a useful tool for nuclear structure calculations of exotic 
configurations, such as halos and extremely elongated fis- 
sioning nuclei. Among the first applications of HFB-AX 
planned are systematic studies of deformed drip-line sys- 
tems. The ability of the code to accommodate very large 
angular momentum cutoffs is crucial in the context of 
nuclei with large neutron skins and halos, in which the 
high-j continuum has significant impact on pairing cor- 
relations [3a |. 

Other applications involve systematic studies of su- 
perdeformed configurations and fission isomers. hfb-AX 
will also provide systematic energy corrections at large 
deformations, essential for HFB models of fission based 
on HO expansion. The differences seen in Fig. [7] are ex- 
pected to appreciably impact the predicted spontaneous 
fission half-lives. 

In this work, reflection-asymmetric and triaxial defor- 
mations, important in realistic fission calculations, have 
not been investigated. We are currently developing a 
symmetry-free coordinate-space 3D HFB solver based on 
the multi-resolution wavelet method. The hfb-AX code 
reported in this paper will provide crucial benchmark 
tests for the general-purpose wavelet HFB framework. 
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